Mixed-precision solver #6521
Mixed-precision solver #6521
Conversation
blattms
added
the
manual:enhancement
|
jenkins build this please |
|
i see mike drop. |
|
1.6 speedup is very promising @nrseman, great work!
|
I have renamed all |
|
Thank you for your interest and your questions, @multitalentloes. I don't want to start a detailed technical discussion in this PR request. I'll contact you separately. |
multitalentloes
left a comment
multitalentloes
left a comment
There was a problem hiding this comment.
Thank you for the contribution to OPM, I have some surface-level feedback on the code right now, and also taken the opportunity to note down some of the goals for this code further down that line.
Immediate code-fixes:
- Compilation error when compiling with CUDA/HIP (see comment)
- Compilation error type mismatch during variable assignment (see comment)
- Ensure that this PR does not break OPM compilation for ARM CPUs (as they do not support avx2)
General stuff:
- Use docstrings, see examples from other files in the codebase
- Use
clang-formatto make it more consistent with the rest of the codebase - Use
fmt::printinstead ofprintf
Future work:
- Registering the preconditioner in the factory (currently trivial placeholder). This will make it available to other multi-stage preconditioners in OPM, for instance letting us combine it with the default CPRW preconditioner, which is would have a great impact if performance is improved! From what we discussed outside of this PR this preconditioner is implemented with knowledge of how the BiCGSTAB is performed and that is something that should be looked more into.
- Having these preconditioners registered in the MPI factory as well with wrapBlockPreconditioner so it can be used in parallel simulations
| DUNE_UNUSED_PARAMETER(prm); | ||
| return std::make_shared<MultithreadDILU<M, V, V>>(op.getmat()); | ||
| }); | ||
| F::addCreator("mixed-ilu0", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) { |
There was a problem hiding this comment.
this line produces compilation warnings due to unused variables
There was a problem hiding this comment.
Addressed by adding
DUNE_UNUSED_PARAMETER(op);and will be included in my next update.
| DUNE_UNUSED_PARAMETER(prm); | ||
| return std::make_shared<TrivialPreconditioner<V,V>>(); | ||
| }); | ||
| F::addCreator("mixed-dilu", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) { |
There was a problem hiding this comment.
unused variable warnings here as well
There was a problem hiding this comment.
Addressed by adding
DUNE_UNUSED_PARAMETER(op);and will be included in my next update.
| virtual void update() override {}; | ||
| virtual bool hasPerfectUpdate() const override {return true;} | ||
| virtual void pre (X& x, Y& y) override {}; | ||
| virtual void post (X& x) override {}; |
There was a problem hiding this comment.
unused variable warnings for all of these functions with arguments as well with the empty body. I think this has been resolved other places in the code with [[maybe_unused]] or a similar attribute
There was a problem hiding this comment.
Resolved by prepending all arguments by [[maybe_unused]] and will be part of my next update.
| printf("n = %lu\n",x.dim()); | ||
| } | ||
|
|
||
| void exportSparsity(const char *path='.') |
There was a problem hiding this comment.
causes compilation error '.' is a char, replace with "." to make it a string/char*
There was a problem hiding this comment.
Interestingly, I don't get a compliation error, but your observation is correct and I have made the fix.
opm/simulators/linalg/mixed/prec.c
Outdated
|
|
||
|
|
||
| /* | ||
| void mat3_view(double const *M, char const *name) |
|
|
||
| } | ||
|
|
||
| int bslv_pbicgstab3m(bslv_memory *mem, bsr_matrix *A, const double *b, double *x) |
There was a problem hiding this comment.
In general for new functions: Use docstring to document what the function does and what its expected inputs, outputs, and special assumptions are to make the implementation more accessible
There was a problem hiding this comment.
@multitalentloes: I'm a bit confused about what you mean by docstring. Do you simply mean a multiiline comment or is it required to adhere to the doxygen standard? I see a mix in the current code. Also, is there a policy on whether the documentation is attached to the function definition of the function declaration?
There was a problem hiding this comment.
Ideally it should adhere to the doxygen standard for those who use that. I would guess declaration is best, at least that is what I have done when declaration is in a separate headerfile
opm/simulators/linalg/mixed/prec.c
Outdated
|
|
||
|
|
||
|
|
||
|
|
There was a problem hiding this comment.
Use clang-format on all the files to ensure that things are compatible with existing code, easiest to just to it now from the start to avoid for instance weird whitespacing such as these blank lines
| { | ||
| public: | ||
|
|
||
| MixedSolver(M A, double tol, int maxiter, bool use_dilu) |
There was a problem hiding this comment.
I am a bit dubious to taking a copy of A here, it seems to work for some inexplicable reasons, but why not take a const ref (const M& A)? I think as it stands, the code is storing a pointer (data_) to a memory location that will be deleted after invocation. No idea why this works in practice here, though.
There was a problem hiding this comment.
The fact that this works implies that a copy of the matrix object A is a shallow copy. Nevertheless, your solution is better and the change will be included in my next update.
|
I have addressed all trivial issues and pushed the updates to the PR branch. What remains is the following:
|
|
@blattms: Would the file below work for the compilation test we discussed? I have verified that the code compiles with the following command gcc -mavx2 -c test_avx2.ctest_avx2.c#include <immintrin.h>
typedef
struct bsr_matrix
{
int nrows;
int ncols;
int nnz;
int b;
int *rowptr;
int *colidx;
double *dbl;
float *flt;
} bsr_matrix;
void bsr_vmspmv3(bsr_matrix *A, const double *x, double *y)
{
int nrows = A->nrows;
int *rowptr=A->rowptr;
int *colidx=A->colidx;
const float *data=A->flt;
const int b=3;
__m256d mm_zeros =_mm256_setzero_pd();
for(int i=0;i<nrows;i++)
{
__m256d vA[3];
for(int k=0;k<3;k++) vA[k] = mm_zeros;
for(int k=rowptr[i];k<rowptr[i+1];k++)
{
const float *AA=data+9*k;
int j = colidx[k];
__m256d vx = _mm256_loadu_pd(x+b*j);
vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+0))*_mm256_permute4x64_pd(vx,0b00000000); //0b01010101
vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+3))*_mm256_permute4x64_pd(vx,0b01010101); //0b01010101
vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+6))*_mm256_permute4x64_pd(vx,0b10101010); //0b01010101
}
// sum over columns
__m256d vy, vz;
vz = vA[0] + vA[1] + vA[2];
double *y_i = y+b*i;
vy = _mm256_loadu_pd(y_i); // optional blend to keep
vz =_mm256_blend_pd(vy,vz,0x7); // 4th element unchanged
_mm256_storeu_pd(y_i,vz);
}
} |
|
@nrseman I think the only thing needed for this to compile on with GPU support enabled is the change in |
|
@kjetilly: After updating my toolchain I ran into the |
Please add a main method to make this compile. We need to be able to run this program. |
There was a problem hiding this comment.
I think there should be some kind of tests that makes sure that the solver runs at least. Currently, we will only know whether it compiles.
If you do C allocations, then you should probably check the results. I am not sure whether the memory is sufficient for arrays that are later realigned.
Is there are any code that checks that:
- the blocks have the correct size
- that the run is serial and not parallel
- etc.
We need to prevent users from misusing this.
I might have missed something of course.
|
|
||
| prec_t *prec_alloc() | ||
| { | ||
| prec_t *P = malloc(sizeof(prec_t)); |
|
|
||
| bsr_matrix* bsr_alloc() | ||
| { | ||
| bsr_matrix *A=malloc(sizeof(bsr_matrix)); |
There was a problem hiding this comment.
the result of malloc should probably be checked.
|
|
||
| bslv_memory *bslv_alloc() | ||
| { | ||
| bslv_memory *mem = malloc(sizeof(bslv_memory)); |
There was a problem hiding this comment.
Might be good enough for OPM, because we never compile with -DNDEBUG to not make our simulator to fast 😅 One day that might change, though.
For "normal" projects I would use a check that even triggers in production. Asserts should be used to catch programming errors. This seems to be something that happens if system resources are exhausted.
| { | ||
| int nrows; | ||
| int ncols; | ||
| int nnz; |
There was a problem hiding this comment.
Does using int for nnz limit us too much?
There was a problem hiding this comment.
Once we start running models with more than 2 billion connections (~25M cells) then yes.
Let us just not lose sight of the fact that we are trying to get a fundamentally new feature into the simulator. Given that the current implementation only allows sequential runs, it is very unlikey that anyone will try to apply it to extrelemy large models.
|
Unfortunately, you will need to rebase this because of conflicts, |
|
This module might get used by others. Please make sure that the new headers are installed by listing them in CMakeLists_files.cmake variable PUBLIC_HEADER_FILES |
|
I am also getting some warnings: I have not checked this and it might be false postives on architecture ppc64el. This is for a rebased version. |
|
It seems like I can call the simulator without --matrix-add-well-contributions=true and do not get an error. What does happen in this case? Maybe we should enforce that option somehow? |
opm/simulators/linalg/mixed/prec.c
Outdated
| _mm256_storeu_pd(xi,vz); | ||
|
|
||
| //double z[4]; | ||
| //_mm256_store_pd(z,vz); |
There was a problem hiding this comment.
some lines of dead code that look like they can be removed
There was a problem hiding this comment.
I tend to leave code snippets in because they are alternative ways of doing the something. It helps me when I go back to the code after a while. Nevertheless, I have removed it just for you. I have removed other commented out code as well.
| for(int k=0;k<9;k++) C[k]=M[k]; | ||
| } | ||
|
|
||
| void mat3_matfms(double *C, const double *A, const double *B) |
There was a problem hiding this comment.
this and other matrix helper functions are not documented
Maybe also extract utility functions as this one to a separate file as it is not really specific to the precs themselves?
There was a problem hiding this comment.
I will not extract the utility functions into a separate file for the following reasons:
- scope creep: This was not a part of your original review and not essential for merging.
- scope: They are currently designed for inline use in the preconditioners, and I do not intend to support general use of them at this point.
There was a problem hiding this comment.
Please document with the mathematical expressions that the functions represent to help other people. Also for the other functions. This will save people a lot of time.
There was a problem hiding this comment.
I have added documentation of to the above function and other local functions.
|
Sorry for the delay. Here is an updated list of things that should be in place before merging:
With these resolved then we can get this preliminary version merged that should be expanded upon in the way discussed earlier in this PR. |
The UMFPackVectorChooser in dune-istl is only specialised for double and std::complex and not float and the template class chosen misses the typenames. We need to work around this here. The question is why the multisegment wells are used with float... |
|
jenkins build this serial rocm hipify please |
|
BTW on my system this PR compiles. I am checking in a chroot without avx support now. |
Maybe this problem is also in master? |
For what it's worth, the CI system builds the master sources with |
Which implies it is not a problem in master. I just compiled my branch against the tip of the current master branch of |
|
Note that the latest version still compiles with |
How did it go? |
Compilation works as expected. |
|
jenkins build this serial rocm hipify please |
Neither I nor jenkins are seeing this. Hence this is not an issue for merging. If your are still seeing this, then please open an issue and provide enough information to replicate |
@multitalentloes This is not even possible on my system: |
opm/simulators/linalg/mixed/prec.c
Outdated
| for(int k=0;k<9;k++) invA[k]=M[k]/detA; | ||
| } | ||
|
|
||
| inline void vec_copy9(double *y, double const *x) |
There was a problem hiding this comment.
Besides all this nitpicking above there is a serious issue. The code does not link because of inline functions when I compile with -O0:
[100%] Linking CXX executable bin/flow
/usr/bin/ld: lib/libopmsimulators.so.2026.04: undefined reference to `vec_copy9'
collect2: error: ld returned 1 exit status
make[3]: *** [CMakeFiles/flow.dir/build.make:233: bin/flow] Fehler 1
make[2]: *** [CMakeFiles/Makefile2:6973: CMakeFiles/flow.dir/all] Fehler 2
make[1]: *** [CMakeFiles/Makefile2:6980: CMakeFiles/flow.dir/rule] Fehler 2
make: *** [Makefile:3195: flow] error 2
| inline void vec_copy9(double *y, double const *x) | |
| void vec_copy9(double *y, double const *x) |
or find another way to include the symbol. The inline is only a hint. Is it really needed? https://stackoverflow.com/questions/19068705/undefined-reference-when-calling-inline-function
The other inline functions here have no problem. That might be because they are actually not inlined?
There was a problem hiding this comment.
Besides all this nitpicking above there is a serious issue. The code does not link because of inline functions when I compile with
-O0
So far, I have not been able to reproduced the issue on my end. Do you mind trying
static inline void vec_copy9(double *y, double const *x)It should limit the scope to the file where the function is defined. I'll keep trying to reproduce the issue.
There was a problem hiding this comment.
@blattms: I've been able to reproduce the -O0 linking issue and verified that static inline fixes it. It will be among the changes I will push shortly.
|
I just pushed an update that take care of your remaining concerns, i.e.
My understanding is that the AVX2 compilation issue taken care of (right Markus?) That leaves us with two items:
So what do you think? Are we ready for merging? |
|
Things are okay from my side now for this version of the code. |
|
jenkins build this please |
blattms
left a comment
blattms
left a comment
There was a problem hiding this comment.
Looks good now. I retested the debugging build. In it goes.
With a light protest because of misusage of assert, though 😅
I'll move away from them going forward. Thanks for all your support! |
7 participants
This PR introduces mixed-precision building blocks for Krylov subspace methods and a highly optimized mixed-precision implementation of ILU0 + BiCGSTAB. Initial testing indicate 1.6x speed-up per linear iteration with minimal impact on iteration count. For information on how to activate the new solver see the README.md file in the
opm/simulators/linalg/mixedfolder.Looking forward to your feedback.
PS! After rebasing on master yesterday, the SPE10 case I used as one of my performance tests failed endpoint validation due to negative SOGCR. I have not investigated further, but it seems to be a finite precision issue. You can see my hack to work around the issue in the
mixed-devbranch on my fork